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We present a theoretical investigation of a lattice Tonks-gas that is created by inelastic, instead of 
elastic interactions. An analytical calculation shows that in the limit of strong two-body losses, the 
dynamics of the system is effectively that of a Tonks gas. We also derive an analytic expression for 
the effective loss rate. We find good agreement between these analytical results and results from a 
rigorous numerical calculation. The Tonks character of the gas is visible both in a reduced effective 
loss rate and in the momentum distribution of the gas. 



I. INTRODUCTION 

A Tonks-Girardeau gas is a one dimensional (ID) gas of 
bosons where a repulsive interaction dominates all other 
energy scales. It was predicted long ago [li|2j that repul- 
sion can mimic the Pauli exclusion principle and that, for 
one dimension, the wavefunction for bosons and fermions 
would be the same, except for a simple transformation. 
Not only that, but the excitation spectrum of these one- 
dimensional hard core bosons can also be mapped to that 
of a fermionic system. A few years ago, experiments ob- 
served a Tonks-Girardeau gas using ultracold gases in 
optical lattices with strong repulsive interactions [3 IH [S] . 

In a recent experiment [6] we showed that strong dis- 
sipation in the form of two-body losses can also simulate 
a Pauli exclusion principle, fermionizing a system and 
transforming it into a dissipative but long-lived strongly 
correlated gas. This equivalence was demonstrated for 
molecules of ^''Rb loaded in an optical lattice. For deep 
enough lattices, these particles exhibit strong two-body 
decay rates which, in appropriate units, exceed the av- 
erage kinetic energy of the particles. The particles then 
avoid coming close together and behave like impenetra- 
ble bosons. This happens both in the continuum case of 
ID tubes and when the tubes are modulated by a deep 
optical lattice — two situations which resemble the ex- 
periments with elastically interacting bosonic atoms in 
Refs. mis] and Ref. [5], respectively. 

At least in the lattice case, the equivalence between 
strong dissipation and a Pauli exclusion principle can be 
understood in terms of the Zeno effect. Following the dis- 
cussion in Ref. [gj , Fig. [l] depicts a stable configuration 
with two molecules on neighboring sites. This config- 
uration is connected via hopping to states with double 
occupancy. These states decay at a rate F which is much 
larger than the hopping amplitude J /h of the particles. 
Treating this as a typical three-level system from Quan- 
tum Optics, one concludes that particles stay on their 
original sites with only a minor loss rate of 0{J^ /h^T). 

In this work we present a rigorous theoretical analysis 
of this system. The outline is as follows. We begin in 
Section |ll] by introducing a master equation that models 
two-body losses for a single species of bosonic particles 
-for instance the molecules in Ref. [Hj-. We will partic- 




FIG. 1: Two molecules sit on neighboring lattice sites. The 
molecules can tunnel with a hopping amplitude J/h which is 
much weaker than the decay rate of the resulting doubly oc- 
cupied site into the vacuum. The effective model is that of im- 
penetrable particles decaying with a much weaker 0{J^ /h?V) 
decay rate. 



ularize the model to the case in which particles are con- 
fined by an optical lattice and explain how dissipation 
becomes the dominant term. In Section UTTl we will show 
that in the limit of strong losses the master equation can 
be replaced by an effective model in which the rapid two- 
body decay has been eliminated. The dominant terms of 
the effective model are identical to the Hamiltonian of 
an elastic Tonks-Girardeau gas. The residual effect of 
losses is a slow perturbation that can be obtained ana- 
lytically. In Section [TVl we compare both results to exact 
numerical simulations of the full master equation. We use 
Matrix Product Density Operators (MPDO) to ver- 
ify that both the density and momentum distributions 
of the particles closely resemble those of a Tonks gas. A 
second signature of the hard-core boson dynamics is a 
slowdown of the inelastic losses, for which we find good 
agreement between numerical and analytical estimates. 
In the last part of our work we offer further details about 
the methods and derivation used to obtain the effective 
Tonks gas models and the slowdown of losses. Finally, 
we summarize this paper in Section [VI| 
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II. FORMALISM 



A. Master equation 



We model the dynamics of the particles using a Marko- 
vian master equation 



(1) 



where H is the Hamiltonian describing the unitary part 
of the evolution and I? is a dissipator associated with 
the losses due to inelastic collisions. In the absence of a 
lattice potential |H] 

H = J d^x^^H,^ + M|3^ J d3xv|/t2vI/2 (2) 
Im(.g3D) 



Vp 



where Hs{Ti) = —h^V^ /2m + Vtrap(x) is the single- 
particle Hamiltonian and ^'(x) is the bosonic field op- 
erator. The strength of the interparticle interactions is 
gsD = 'inh^a/m, where m is the mass of a particle and o 
is the scattering length. Re(a) describes elastic collisions, 
whereas Im(a) < describes inelastic collisions that lead 
to losses. 

To assert the consistency of this model, let us estimate 
the decay rate of the particles. When taking expectation 
values over the density operator, n{x) = ^''^ (x)^'(x), the 
previous master equation becomes 

|(v|/t(x)vl/(x)) = ^Im(.g3D)(*t'(x)*2(x)). (4) 

For a condensate of molecules this gets the expected form 
d 



dt 



n(x) ~ -i4:3Dn(x) 



(5) 



where n(x) is the density of bosons at a given point and 
the decay rate is proportional to the square of the density 
and to K^Y) — — 2 Im(g3£i)/?i. 



B. Optical lattice 

As explained in the introduction, in this paper we want 
to model the experiment with molecules in a deep optical 
lattice [6 . For low temperatures and tight confinements, 
we can expect that the particles will accommodate to 
the motional ground state of each lattice site. Note that 
if this is true for individual atoms, it is even more so 
for molecules, because having twice the mass and twice 
the polarizability, they are better trapped by the same 
optical potential. Under these conditions, it was shown 
in Ref. (9j that it is convenient to expand the bosonic 



field operator using Fock operators, a^, that create 
particles on the fc-th lattice site: 



^4w(x-Xfc). 

k 



(6) 



The w(x) are Wannier wavefunctions associated to the 
states localized on each site. In particular, for the ar- 
rays of one-dimensional lattices used in Ref. [6] , we can 
separate the wavefunction of these bosonic modes into a 
product of three wavepackets, w^x) = w^(x)w\^i^y)w\_^{z)^ 
a less confined longitudinal one, lyy, and two transverse 
ones, which are very tight due to the perfect decou- 
pling between adjacent tubes. 

The advantage of the Wannier functions is that we 
can now perform a tight-binding approximation: any- 
where outside the single-particle terms in the Hamilto- 
nian, the overlap of different Wannier functions is ne- 
glected. This procedure transforms the unitary part of 
the master equation into a Bose-Hubbard model [9J , and 
discretizes the dissipative terms as well 



{ki) " k 

Y^/r, 2 t2 t2 2 t2 1\ 



(7) 



Vp 



Regarding the notation, the sum (fc, I) extends over near- 
est neighbors along the same tube, |A: — ?| = 1. The tun- 
neling amplitude between neighboring lattice sites is de- 
noted by J. Finally, the on-site interaction matrix ele- 
ment contains both a real and imaginary part 



3zB \ dx\w\\{x)\'' 



dv\wy_{y)\^ 



(8) 



[/,. + iUi ^ Ur — i- 



hT 



which contribute to the unitary and the dissipative parts 
of the master equation, respectively. 

The imaginary part of the interaction constant governs 
the decay of the number of particles per site, rik — a\,ak 



dt 



(uk) = -T(nk{nk ~ 1)). 



(9) 



In the cases that we will study this rate will be larger 
than the speed at which particles tunnel to neighbor- 
ing sites, r ^ J/h. To facilitate the calculations, we will 
group the terms in the master equation according to their 
strength. We introduce a superoperator V that contains 
the tunneling part {Hj cc J) and a superoperator £int 
that describes the elastic {Hc\ oc Ur) and inelastic inter- 
actions {V oc r) 



dt 
Vp 



P = 



(V + /lint) P 

-;^iHj,p] 



:Vp. 



(10a) 
(10b) 
(10c) 
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The term "superoperator" refers to the fact that V, V, 
and Cint are linear operators acting on density matrices, 
not on pure states. Finally, in the following it will be 
important to realize that V is of order J, whereas Cint 
is of order \U\ ^ |J|, and dominates the evolution. The 
following sections exploit this difference of scales to cre- 
ate new and simpler effective models that describe the 
dynamics of the molecules. 



III. EFFECTIVE MODELS 

A. Second order effective theory 

Our goal is to develop an effective master equation 
that is equivalent to (10 1 in the limit of strong dissi- 
pation, TiT ^ J, in which £int dominates. We will 
sketch the main ideas and refer the reader to Sect. |V] 
for more details. The process begins by identifying the 
eigenvalues, A^, and eigenspaces of the dominant term, 
'Cint = then decompose the density matrix 

into a sum of contributions from these eigenspaces 



(11) 



Finally, we will determine some approximate evolution 
equations for the different pi in the presence of a non- 
zero hopping term, V. 

We only need to study three eigenspaces, correspond- 
ing to the dissipation-less states, Aq = 0, and to the states 
most immediately connected to them by the hopping. 
The most relevant term of our density matrix, po 0{l), 
is given by states with zero or one particle per site. These 
states do not decay in the absence of hopping and have 
the property that they belong to a space with zero or one 
particle per site. We thus have 



Po = QopQo 



(12) 



with a projector expressed in the Fock basis of occupation 
numbers 



)o = (|0)(0| + |l)(ir^ = g« 



(13) 



where L is the number of sites in the lattice. The next 
states that we need to consider have a pair of particles on 
some site. As we will see later, we have two sets of states 
depending whether the double occupation is on the left 
or on the right side of the matrix coherences 



VlaP 

VibP 



QipQo, 
QopQi, 



(14) 



where we have introduced 

L 



L-k 



(15) 



In the absence of hopping these states decay with eigen- 
values Ala = -T/2 - iUr/fi = -iU/Ti and Aib = -T/2 + 
iUr/Ti — ill* /U, respectively 

By neglecting higher order contributions to the density 
matrix, it is possible to integrate formally the projected 
master equation and obtain an effective model for the 
dominant term, po{t). After some manipulations one ar- 
rives to the following model [Sect. V C| 



dpo 
dt 
Ci 

Co 



{Ci + C2)po (16a) 

VoWo (16b) 

J2 ^T^oWcWo, (16c) 

ce{la,lb} ^ 



We will now write down explicitly and discuss the mean- 
ing and implications of these different terms. 



B. Tonks-Girardeau Gas 



The most important term in our effective model (16 1 
is given by the Liouville operator Ci = VqVPo- As 



shown in Sect. VD this superoperator is equivalent to a 
Hamiltonian of hard-core bosons, also known as a Tonks- 
Girardeau gas. Therefore, to lowest order 



dt 



Po = 



(kA) 



+ 0{J'/\U\), (17) 



with bosonic hard-core operators cj. and Ck that have 
implicitely the restriction of one particle per site 



Ck - |0)fc(l| 



4 = |l)fc(0|fe. 



(18) 



This is the main result of our paper. Namely, that a 
strong dissipation such as the two-body losses from our 
system can lead to coherent evolution. As it was men- 
tioned in the introduction, the same result can be ob- 
tained in an alternative way. By establishing an analogy 
between losses and a continuous measurement, it is intu- 
itively clear that the strong dissipation causes a Zeno ef- 
fect which suppresses the process of two particles coming 
together and being lost. In this regime, the dynamics of 
the molecules must be given by a Hubbard model where 
doubly occupied states have been projected out. 

In practice, the hard-core bosons model implies a very 
simple dynamics that can be tested numerically, as we 



fc=0 



do in Section IV However, verifying the same thing with 
real particles in an optical lattice represents a challeng- 
ing experiment. Nevertheless, one can easily check two 
phenomena: first, that in the regime of strong decay, 
?iT 3> J, no significant losses take place in a time scale 
of order l/F; and second, that the effective loss rate can 
be estimated using the second part of our effective model 
(16"). This is the goal of the following sections. 
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C. Second order corrections 

With a lengthy calculation we can rewrite the second 
order terms in Eq. ( 16 1 as an effective Liouvillian £2 with 



both a Hamiltonian and dissipation 



i 1 



(19) 



It is convenient to introduce an operator that destroys a 
pair of particles in neighboring sites 

Cfc = Cfc(cfe+i +Cfc_i). (20) 

With these pairs the Hamiltonian part can be written 



k 



4/ 



with an effective strength (1/Aia = hi/U) 

1 \ 2J2 



2.P 
J2 = -^Im 
n 



Ala 



;Ur. 



(21) 



(22) 



Note that this Hamiltonian contains both effective near- 
est neighbor interactions and a three-site hopping of the 
form cj^^^rikCk-i- Both terms are typical of the Bose- 
Hubbard model in the limit of strong repulsive interac- 
tion Ur ^ \J\ [IQj. Moreover, these Hamiltonian correc- 
tions disappear when the lossy particles do not interact 
elastically on-site, Ur — 0. 

The dissipative term is equally simple, 

V2P0 - J2 (^CkPaCl - ClCkPo - PaClCk) , (23) 



and has a loss coefficient 
2J2 

T2 = -^Rc 



1 

Ala 



2.P 



(24) 



In the limit of weak elastic interaction between molecules 
we may write 



2J2 



1 



ul 
u? 



4J2 



1 



ml 



, (25) 



which shows that the decay rate is indeed 0{J^ /Ti^T) as 
anticipated. 



D. Effective losses 

Let us write the evolution of the total number of par- 
ticles, iV, under the effective master equation (16) 



dt 



(N) ^-J2mclCk,N]) ^ -4T2j2{ClCk), (26) 



This expression in general cannot be simplified any fur- 
ther, at least not without some assumption about the 
state with which we compute the expectation value. 



The experiments described in Ref. [B] start from a state 
with exactly one molecule at each site and evolve only 
until approximately half of the particles are lost. Here 
further approximations are possible: first we treat the 
system as homogeneous and second we assume that the 
populations of different sites are uncorrelated[T6]. We 
thus obtain 

(ClCk)^ E {ci4ckci,)^zn', (27) 

l,l'e{k-l,k+l} 

where n is the density and z = 2 is the coordination 
number of our lattice. This approximation leads to a 
rate equation characteristic of two-body processes 



di 



n ~ -4zr2n^ 



-KTl 



(28) 



A similar equation was derived previously in Ref. 6 us- 
ing a more restricted theory than our effective master 
equation ( 16 ). 



IV. NUMERICAL RESULTS 

In order to study the quality of the approximations 
used in the effective model, we have performed numer- 
ical simulations of the full master equation ( 10 1 using 
MPDO 13 • This is a method that approximates the den- 
sity matrix p(t) using a matrix product state structure. 
As described in Ref. [7J, this variational ansatz is well 
suited to simulating evolution of the state under a Liou- 
villian like V -l- >Cint J which can be decomposed into a sum 
of local or nearest-neighbor terms. In our simulations we 
have worked with up to 48 sites and open boundary con- 
ditions, setting [/^ = so that the different effects cannot 
be attributed to a repulsive interaction. We have exper- 
imented with different cutoffs, from two to four particles 
per site, verifying that they give similar results. 



A. Comparison with the Lossless Tonks Gas 

Fig. |2] shows numerical results for TiT / J = 4000. The 
system was initially prepared such that the A^ = 24 cen- 
tral sites of the lattice contain one particle per lattice 
site, while all other lattice sites are empty, similar to the 
state experimentally prepared in Ref. [TT]. We then let 
the system evolve for a time t = 2h/J and measure both 
the density and the momentum distributions, which are 
plotted in Figs. [2](a) and |2jb), respectively. The solid 



lines were obtained with the full master equation (10 1 



and differ only slightly from the dashed lines obtained 
with the lossless Tonks model of Eq. (17 1. The observed 



difference in the distributions is largely due to the re- 
duced particle number. 

In order to quantify how small the remaining difference 
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FIG. 2: (Color online) Comparison between the dissipative 
system (solid lines) and the lossless Tonks gas (dashed lines). 
A'^ = 24 central sites of the lattice initially contain exactly one 
particle per site [dotted line in (a)]. The system is evolved 
with hT/J = 4000 for a time t = 2h/J. The evol ution is sim- 
ulated either with the full master equation ( 10 1 or with the 
lossless Tonks model Eq. ( |17[ |. Part (a) shows the position 
distribution versus site index, part (b) the momentum dis- 
tribution versus quasi-momentum p (a is the lattice spacing 
and L the lattice length), both after the time evolution. The 
difference between the distributions is marginal, except for an 
overall reduction of particle number. 



of the position distributions is, we define 



E 



Tonks 



N{t) N{0) 



(29) 



where N(t) is the total number of particles at a given 
time. Tlx is the number of particles for the site x and 
^Tonks jg same for the lossless Tonks gas. A similar 
measure can be defined in momentum space. Fig. [sj^a) 
shows these quantities. For J/T the distributions 
converge to those of a lossless Tonks g expected 
from the effective model. 

Another important quantity in this comparison is the 
fraction of the state in the hard-core bosons subspace of 
states with at most single occupation 



F = Tr(Qop) 



(30) 



This quantity is shown in Fig. ^h). Comparison with 
part (a) shows that 1 — F tends to be much smaller than 
e. This is because loss causes the system to evolve into an 
incoherent mixture of systems with different total parti- 
cle number. Each of these systems remains to a pretty 




p 0.98 



0.0125 0.025 

j/nr 



0.0375 



FIG. 3: (Color online) (a) Difference between the distribu- 
tions of the dissipative system and the lossless Tonks gas. 
Calculations as in Fig. |2] for a variety of values J/hV. Each 
yields a different distribution for the dissipative system and 
the lossless Tonks gas. The value e defined in Eq. (29 1 quan- 



tifies this difference for the position distribution (solid line) 
and the momentum distribution (dashed line), (b) Fraction 
of the density matrix in the hard-core bosons subspace, with 
F defined in Eq. 1301). 1 — _F is typically smaller than e. 



good approximation in the Tonks gas subspace, but the 
position and momentum distributions look different so 
that the mere scaling with total particle number yields 
larger deviations in e. 



B. Effective Loss Rate 

Our analytic model suggests that the loss rate of par- 
ticle number in a homogeneous system can be approxi- 
mated by a two-body decay Eq. (28 1, which corresponds 



to an evolution of the density given by 



n{t) 



1 + n(0)Ki' 



(31) 



This result applies to the case in which different sites are 
not very much correlated and the losses F dominate over 
the hopping, J. 

In order to test this prediction, we have performed 
MPDO simulation of Eq. ( 10 1 using a uniformly filled 
lattice with L — AQ sites and = 40 particles. We 
studied the evolution for a wide range of values of J/TiT, 
with some examples shown in Fig. [4|a). After an initial 
transient that vanishes on a timescale w l/F, a two-body 
decay Eq. (28 1 fits the data fairly well. However, while 
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FIG. 4: (Color online) We study the density as a function 
of time, for a simulation of Eq. ( 10 1 using L = 40 sites, each 



one initially filled with exactly one particle, (a) Evolution of 
the density as a function of the adimensionalized time tV, for 
three values of the losses J/hF = 0.0179, 0.0541 and 0.1 (from 
top to bottom). We present data from both the numerical 
simulation and for a fit with Eq. ( |34[ |, using circles and lines, 
respectively, (b) Best-fit values of the effective decay rate k 
vs. the value of J/hF at which the simulation took place. The 
solid line represents the theoretical prediction Eq. (281. The 



crosses and the circles corresponds to fits using the formulas 
in Eqs. (|32l) and (|34l), respectively. 



for J/fiT <C 1 the effect of the transient is negligible, 
for larger J/TiT there are better ways to fit the resulting 
curves. One possibility, used in Ref. [B| to fit the numer- 
ical data, is to include a free parameter to describing an 
offset on the time axis 



nit) 



m 



l + n{0)K{t-to)' 



(32) 



However, we found that a more accurate model is a mod- 
ified decay equation with an exponentially modulated de- 
cay coefficient 



dt 



-k(1 — e 



(33) 



Here the exponential term with A F represents 
an heuristic approximation to the transients that we 
have neglected in developing the effective model [See 
Sect. VC . The solution of this differential equation still 



has two fit parameters 



n{t) 



1 + n{Q)K{t + [exp(-At) - 1]/A} ' 



(34) 



Fig. |4]^b) shows the loss rate coefficient k from these two 
fits and Fig. |4|a) shows the quality of the fits for small 
values of the hopping. For larger values, though, the 
fitting becomes numerically unstable. In this regime a fit 
(32 1 behaves better. 



What we do not show in the previous plots is that the 
long term behavior of the system no longer follows the 
simple two-body decay laws from Eqs. ( 32 1 and ( 34 1 . The 



reason is that at low densities there are enough correla- 
tions that we can no longer use the simple models from 
Section [III C[ In this regime ofhe small densities, a coarse 
grained description becomes approximately equivalent to 
the Lieb-Liniger model [12j but with inelastic interac- 
tions. As we have shown elsewhere [5], the decay at long 
times is then expected to follow the law dn/dt cx h'^. 



V. DETAILED CALCULATIONS 

The goal of this work is to find an effective model for 
the particles in the lattice, which works in the limit of 
fast dissipation, J <C KT. Our main tool to understand 
this limit is a generalization of Kato perturbation theory 
|13] . also known as adiabatic elimination, to the superop- 
erators V and £int- Section |VX] explains how our calcula- 
tions relate to this broader scope. Readers less interested 



in this aspect may skip to Section VB where the actual 
derivation begins. 



A. General ideas 

Both V and Ci^t Siie linear operators that act on an 
appropriate space of matrices. Even though within this 
space the superoperators are not Hermitian, the domi- 
nant term £int has infinitely many eigenstates forming a 
discrete spectrum of well separated points beginning at 
Aq = and spanning through the left half of the com- 
plex plane [Fig. [5|d]. Each of these points represents a 
family of matrices that, in the absence of hopping term, 
decay at a rate given by the real part of the eigenvalue, 
Re(Ai) < 0. In the limit of strong dissipation, the hop- 
ping V will couple these eigenspaces very weakly, with an 
amplitude J/h much smaller than the typical eigenvalue 
separation, F. 

Precisely in this limit of weak hopping we will be able 
to write analytic expansions of the eigenstates, eigenval- 
ues and evolution equations for the perturbed superop- 
erator (V -I- Cint)- The ideas is to use Kato's resolvent 
method |13J with an expansion of the Liouvillian 



•^int — AiT^i, 



(35) 



which uses a complete set of pseudo-projector operators 
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FIG. 5: (a) In ordinary perturbation theory of Hermi- 
tian operators we find eigenspaces with eigenvalues which 
are well separated along the real axis. A weak Hermitian 
coupling, smaller than the energy separation between spaces 
e ^ |Ai — Aj|, only causes small shifts in these energy levels, 
(b) Our Liouville operator C\nt is not Hermitian but has well 
separated eigenvalues in the space of density matrices (We 
have set Ur = G for simplicity). The real part of these eigen- 
values now represents the decay rate of such matrices. The 
hopping of particles J acts as a weak perturbation J ^ KV, 
that slightly changes both the real and the imaginary part 
of the original eigenvalues. In particular, the states with one 
particle per site, po, now acquire contributions of 0{J /VY and 
(!?( J/F)^'^' " by coupling directly, pi, or indirectly, p2,3,..., to 
faster decaying subspaces. Both in (a) and (b), for second 
order approximations of the unperturbed eigenspaces, Ao, we 
may neglect all spaces with an indirect coupling. 



These operators, built from the right and left eigenvec- 
tors of >Cint, are non-Hermitian but this fact poses no 
difficulties in generalizing the usual perturbation theory. 

In particular, the space of density matrices that had 
zero decay rate for J — Q will transform into the per- 
turbed eigenstates and eigenvalues 

(V + Ant)po( J) = Ao(J)po(J). (37) 
The decay rate will no longer be zero but remain small 

J 



Ao(J) c„ 



n>l 



o{r/nT), (38) 



and these states will acquire a small admixture of the 
matrices belonging to the unperturbed eigenspaces that 



had larger decay rates 



n>l ^ '' 



(39) 



Notice the weight of the different contributions pn de- 
pends on how many applications of V we have to perform 
to connect po to these other eigenspaces. Furthermore, 
we understand that it is the coupling of po to the de- 
caying states through V what makes the initially stable 
space lose particles. A similar reasoning can be applied 
to the rapidly decaying eigenspaces, pi, 2,. ..(•/), which will 
be modified by the coupling V. However in this case the 
decay rates will remain of order F with small corrections 
from the hopping of particles. 

We therefore conclude that if we prepare our particles 
in any initial state, after a short transient t 1/r, most 
of the state will be captured by the eigenspace with the 
lowest decay rate (38). Following the structure given in 
Eq. ( 39 1 we will decompose our evolved state in a series 



of contributions from the unperturbed eigenspaces 

p{t) ^ P0{t) + pi{t) + p2{t) + . . . , 



(40) 



where, according to Eq. ( |39[ ), the high order contributions 

have a vanishingly small size 



|p„(t)| ~ (JAF)"|po(i)l, t->J. 



(41) 



We are therefore allowed to perform a self-consistent ap- 
proximation which consists on, first, writing the evolu- 
tion equations for each subspace Pnit), then impose that 
all contributions Pn>2 — and finally integrate the re- 
maining equations until we reach an effective model for 
the lowest order term. 



B. Local Projections 

As sketched before, our study of the losses in the lattice 
begins by finding out the eigenspaces of the unperturbed 
superoperator, £int. This task is greatly simplified by the 
fact that Cint — X^/f^ioc.fc is a sum of commuting local 
superoperators >Cioc- We can thus focus on diagonalizing 
one of these superoperators on a single lattice site. 

We will now introduce some notation. Since we are 
interested in density matrices as elements of a Hilbert 
space on which the superoperators act, we will introduce 
a basis of vectors of this space. For a single site we will 
define the Fock states 



|n, m) 



1 



(42) 



and introduce the Hermitian conjugate (n',rn'| given by 
(n' ,m'\n,m) = S„n'Smm'- In this basis C\oc becomes a 
bidiagonal, nonsymmetric operator 

-^(Cn + U)\n, m) + ^ ye^ln - 2, m - 2),(43) 
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where = x{x — 1). The kernel of this operator is given 
by states with or one particle per site, so that we can 
write the space of non-decaying density matrices as in 
Eq. pi). 



In addition to the right eigenvectors, (>Cioc ~ A„)|u„) — 
0, we will also search the left eigenvectors, (w„|(£ioc — 
A„) — 0, which have the property {wn\vm) = Snm- With 
these families of operators we will construct a set of 
pseudo-projector operators of the form 



tIoc 



\Vn){Wn 



(44) 



We face one problem, though, which is that £ioc acts 
on an infinite-dimensional space, with occupation num- 
bers that can be arbitrarily large. We argue that for our 
purposes it suffices to truncate the Hilbert space to occu- 
pations n,m < 2. The reason is that our initial states will 
all belong to the space po given above and, since we will 
neglect contributions from third and subsequent order 
couplings [i. e. pi.z,... = in Eq. (40 1], we will obtain at 



most a double occupation per site. The self-consistency 
of our approximation will be evident at the end. 

Using the previous truncation, C\qc becomes a 9 x 9 
block-diagonal and banded matrix with projections onto 
eigenspaces 



^0 



ploc 
' la 



-ploc 
' lb 



&,f)'=0,l 



|0,0)(2,2H- M{b,h'\ 

b,b' 

6=04 

E l^'2)(^2| 



6=04 

|2,2)(2,2|-|0,0)(2,2| 



and corresponding eigenvalues 



Ac 


= 




Ala 


iU 


r 


~T ^ 


^2 


-^16 






A2 


= -r. 





iUr 
IT 



(45a) 
(45b) 
(45c) 
(45d) 

(46a) 
(46b) 

(46c) 
(46d) 



and Vl°i,'^ are connected to Pq°'^ by a single appHca- 
tion of the hopping term, V, and 7^2°'^ by two applications 
of V. Note also that V^"^ contains not only a projector 
onto single or zero occupancy, but also a term that de- 
stroys two particles at a site, emptying it. This will be 
essential later on. 



Adiabatic elimination 



As explained in Section III we will consider that our 
states can be described by a density matrix with con- 
tributions coming from the eigenspace with the lowest 
decay rate, po, and the two eigenspaces connected to it. 



Pia.ib- These contributions are respectively obtained by 
applying the following pseudo-projectors onto p{t) 



-Pla - E(^^) 



L-1 



loc\ {g>L— m- 



(47a) 
\(47b) 



Pib 



m=0 
L-1 

E 

m=0 



loc\ ®m 



'^1°6'® (^o°')^^"'""'-(47c) 



Connecting to the previous notation of projectors onto 
states with zero and one particles per site, it will be useful 
to realize that the zeroth order projector can be written 
as follows 



'Pop = QopQo + I E '^iQiPQi^ 



,t2 
'k ■ 



(48) 



Since we will neglect higher order couplings, we can 
use Eqs. (10a I and (36 1 to write evolution equations for 



the density matrices in the form 

dpo 
dt 



VooPo + E 



dpc 
dt 



= XcPc + VcOPO- 



(49a) 



(49b) 



We have abbreviated Vij = ViVVj and introduced the 
notation that the index c runs through {la, 16}. The 
terms Pc can be integrated out of the model using the fact 
that our states are initially prepared in the slow decay- 
ing manifold Pc{0) ~ 0. Formal integration of Eq. (49b I 
yields 



Pc{t) 



dre 



VcoPo(t), 



(50) 



which after integration by parts becomes 



Pc{t) 



t'-Vco [po{t) 

Ac 



A. 



Po(0)] 
dpo 



dre-^=-Vco^(r). (51) 



dt 



We neglect the remaining integral, because it is of higher 
order in J/fiT than the previous term, which is evident 
from the fact that dpo /dt oc J in Eq. ( 49a I . Insertion of 
Pc{t) into Eq. (49a I yields 



dpo 
dt 



- Voo- 



J2 Y^OcVco] Poit) 



y- 



VocVcO/5o(0). 



(52) 



The first line of this equation represents our effective 
model ( 16 1 and will be discussed in the following sec- 



tion. The second line is a transient that decays at a rate 
cx r. Therefore we obtained that the system converges 
to the slow-decaying eigenspace in a time t 1/r, much 
shorter than the typical time scale fi/J at which the ef- 
fective model operates. 
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D. Hard-core bosons 

Let us analyze the lowest order contribution to our 
effective model ([52| , given by £i = Vqo- Using the ex- 
pression in Eq. ( |48[ ), we obtain 



Qo^[Hj, QqPqQo]Qo 
-^[QoHjQo,po]. 



(53) 



fn other words, this Liouville operator is equivalent to 
a Hamiltonian in which we have projected out all states 
with double occupation. This is the hard-core bosons or 
Tonks-Girardeau gas model presented in Sect. |fff B| 



E. Second order losses 

We are now going to consider the second order Liou- 
villian £2 from Eq. ( 16 ) 



E 



-1 



A, 



(54) 



ce{la,lh} 

We can expand this expression 



It is now time to rewrite everything in terms of hard core 
boson operators. We notice the following equivalence 



-,t2 



{k.l) 



k 



(58) 



This arises from the fact that Qi projects onto a state 
with a single pair. Therefore, the tunneling term only 
contributes with processes that take two neighboring par- 
ticles (Cfe) and create a pair in one of the sites. Notice 
that Cfe already enforces the projection Qq- Using this 
notation we can simplify our expressions even further 



4t 



(59) 
(60) 



thus arriving to the final model 



C2P0 - ^2 



2J2 



2 2Re-^^C,poCl, (61) 

Ala , 



which is studied in Sect. IIII CI 



C2P0 = 



+ 



Xibh' 



:ro[Hj,Qi[Hj,po]Qo] 



(55) 



:Vo[Hj,Qo[Hj,po]Qi]. 



Using the property po = QaPaQa, one realizes that the 
only relevant terms are £2^0 = 'Po-^Po/f^^ with 

1 



Apo 



1 



lb 



{HjQiHjQoPo — QiHjQqPqQoHj) 
{QoPoQoHjQiHj - HjQoPoQoHjQi^^e) 



We now consider the final projection with Vq. Following 
Eq. (48 1, this pseudoprojector contains two operations: 



the first one keeps terms proportional to QqHjQiHjQq, 
while the second one acts on the terms that create a dou- 
bly occupied site on each side of the density matrix, that 
is QiHjQopoQoHjQi. Introducing T = QiHjQq/{-J), 



^2Po = ^ i^T^Tpo + —po 
2J2 



1 TtTl 



(57) 



VI. CONCLUSION 

We have shown analytically and confirmed numerically 
that strong, inelastic interactions can induce a Tonks gas 
dynamics for a cloud of molecules trapped in an opti- 
cal lattice. The particles act like hard-core bosons, with 
dissipation playing the role of a strong repulsion. This 
effective model is completed with a reduced loss rate, 
7eff oc J^/r which is much slower than both the tunnel- 
ing amplitude, J/h, and the original loss rate, T. 

Even with the small losses, the state of the system 
can at all times be described as an incoherent mixture of 
strongly correlated Tonks gases with different total par- 
ticle number. In this respect, being based on the idea of 
using dissipation to create strong correlations, our paper 
connects to recent works which suggest using dissipation 
to engineer states and phase transitions [l4l [TS]. 
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